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ABSTRACT 

In this paper, we present a non-oscillatory spectral Fourier method for the solution 
of hyberbolic partial differential equations. The method is based on adding a nonsmooth 
function to the trigonometric polynomials which are the usual basis functions for the 
Fourier method. The high accuracy away from the shock is enhanced by using filters. 
Numerical results confirm that no oscillations develop in the solution. Also, the accuracy 
of the spectral solution of the inviscid Burgers equation is shown to be higher than a fixed 
order. 
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1. INTRODUCTION 


In this paper, we discuss shock capturing techniques using spectral methods. In par- 
ticular, we would like to present a nonoscillatory version of the spectral Fourier method 
when applied to a nonlinear hyperbolic equation. The main difficulty in applying spec- 
tral methods to discontinuous problems is of course the Gibbs phenomenon. In fact, this 
problem exists even in the approximation level. It is well-known that if a discontinuous 
function f(x) is approximated by its finite Fourier series P^f 

P K f = £ ?*«“•, 

k=—N 

/fc = ^ Jo f{x)e~ xkx dx , (1.16) 

then the order of convergence of P^f to / is only O(-jk) for each fixed point. Moreover, 
Puf has oscillations of order 1 in a neighborhood of O(-j^) of the discontinuity. 

In the applications, we usually have piecewise C°° functions, and in this paper we will 
consider only those functions. It is known that it is possible to improve the accuracy of the 
approximation away from the shocks. There are currently two methods (see [7]) that are 
being used. The first amounts to modifying the Fourier coefficients by multiplying them 
by a decreasing function cr(A;) . Some of the commonly used filters are 

Ok = e -a ^~7r IL ^ ra \k\ > ko 

( 1 . 2 ) 

o k = 1 |fc| < k 0 . 

The second method [l] is based on convoluting the approximation with an appropriate C°° 
function ip(x,y) such that 


Pxf • >/>(i,y) ~ /(</)■ 


(1.3) 


While both (1.2) and (1.3) are effective away from the discontinuity, they do not eliminate 
the Gibbs phenomenon in the neighborhood of the shock. This is very important for the 
stability of the spectral method when applied to partial differential equations. In fact in 
Section 2 we show that the total variation of P^f grows like logiV. It is easily shown that 
this is the case also for the filters in (1.2). 
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In Section 2, we show that by adding a sawtooth function to the basis functions e' kx 
one can control the Gibbs phenomenon. This in conjunction with the filters (1.2)-(1.3) 
yields a high order nonoscillatory approximation to a piecewise C°° function. In Theorem 
2, we prove that the total variation of the new approximation converges to that of the 
approximated function. We also prove that the convergence for the new approximation in 
L 1 norm is one order higher. 

Many modern nonlinear schemes for the solution of the conservation equation 

«* + /(«)* = 0 (1.4) 


are based on two distinct steps, namely reconstruction and time marching. We use the cell 
averaging formulation to rewrite (1.4) as 


duj 1 
dt Ax ; - 


{f ( u .+i) - f (“»-d) 


= 0 


(1.5) 


where 


A x i ~ x i+ k x i~ h' u *~ As, L ** 

3 i~* 


3+1 udx 

h 


f{u j+ i) = f{u{x J+ l)). 

The first step, then, is to reconstruct the function it(s) from u(s) . It is here that we 
use the nonoscillatory technique developed in Section 2. For the second step, the time 
marching, we use the third order Runge-Kutta scheme developed in [12]. We try to avoid 
any modification technique (like the application of limiters) in order to avoid deterioration 
of the overall accuracy. 

We demonstrate in the last section that the procedure applied to several model problems 
yields indeed nonoscillatory results with an order of accuracy which is higher than algebraic 
away from the discontinuity. 


2. NON-OS CILLATORY APPROXIMATION 

In this section, we suggest a method to reconstruct a nonoscillatory approximation to a 
piecewise C°° periodic function from its first N Fourier coefficients. The approximation is 
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nonoscillatory in the sense that the total variation of the approximation converges to the 
total variation of the approximated function. Moreover, the approximation convergences 
in the maximum norm outside a small interval around the point of discontinuity. Applying 
the filters (1.2)-(1.3) will increase the order of convergence away from the discontinuity, 
thus providing a nonoscillatory spectral approximation. 

For simplicity, assume that u(z), 0 < x < 2tt is a periodic piecewise C°° function with 
only one point of discontinuity at x, and denote by [u] the value of the jump of u[x) at x s , 
namely 


r , _ u(x+)-u(o;J 
2tt 

We assume also that the first 2N+ 1 Fourier coefficients ut of u(x) are known 


( 2 . 1 ) 


u t = 


2tt 



u(x)e ,tx dx, 


-N<l<N. 


( 2 . 2 ) 


The objective is to construct a non-oscillatory spectrally accurate approximation to u(x) 
from the U( s. We start by noting that the Fourier coefficients u*’s contain information 
about the shock position x s and the magnitude [u] of the shock. In fact we can state 


Lemma 1: Letu(x) be a periodic piecewise C°° function with one point of discontinuity 
x si then for \l\ > 1 and for any n > 0 


V' 1 + J_ F 

t' 0 (H) k+1 Jo [MY 


u t = e 


(2.3) 


Proof : Since 

1 f 2t •» 1 r x ‘ 1 r 2r 

tit = — / u[x)e ltx dx= — / u[x)e ,tx dx+ — / u[x)e ,tx dx 

27 r Jo 27 r Jo 27r Jx , 

we can integrate by parts to get 

u K)~ u ( x 7) , J_ [ 2 * u '( x ) e ~ itx dx 
Ut ~ 6 2ml ' 27 t Jo M dX 

the rest is obtained by induction. This completes the proof. 


(2.4) 
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As an example, consider the sawtooth function F(x,x s ,A) defined by 

F(x,x„A) = A | 2 *_ 


x < x. 


(2.5) 


X X> X, 

Note that the jump of the function - [F] - is A and all the derivatives are continuous 
[^r] = 0 for k > 1. That means that the expansion (2.3) can be terminated after the first 
term, yielding the following results for f k - the Fourier coefficients of F(x,x a ,A) 


-ikx. 


/*(„,, A)= A— 


1*1 >1 


ik 

fo{x 3 ,A) = A{tt - x a ). 

This example suggests that we can rewrite (2.3) as 

ft, = /<(*., H) + E ^ C * ix ' |£| - L 

The order 1 oscillations in approximating u{x) by its finite Fourier sum 


( 2 . 6 ) 


(2.7) 


N 


Pnu = Y u t e' 

t=—N 


tlx 


are caused by the slow convergence of 


N 


F N (x,x„ [tz]) = Y Mx„[u})e llx 
t=-N 


( 2 . 8 ) 


(2.9) 


to the sawtooth function F(x,x a ,[u]). Therefore, those oscillations can be eliminated by 
adding a sawtooth function to the basis of the space to which u(x) is projected. To be 
specific, we seek an expansion of the form 


v N (x) = Y <h* itx + Y -r p e-^e itx 

|*|<AT \1\>N 


( 2 . 10 ) 


to approximate u(x). The 2 N + 3 unknowns a*(|£| < N), A and y are determined by the 
orthogonality condition 

ri* 


r 2*- 

/ (u — Vtf)e~ ,kx dx = 0 |fcj < iV + 2. 
Jo 

The system of equations (2.11) leads to the following conditions: 


( 2 . 11 ) 


a, = ft, [£| < N 


( 2 . 12 ) 
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(where u t are the usual Fourier coefficients of u(x), see (2.2)) and 


i[N + 1) 


<r*' (W+1)y = u N+ 1 


A c ~ i(N+2) V _ A 

i(N + 2) 6 ~ UN+2 ' 


Solving (2.13) for A and y one gets 


e' v = 


_ (-W + l)u^+i 
(N + 2)ujv + 2 


\A\ — (N + 1)|ujv + i|. 


(2.13a) 

(2.136) 


(2.14a) 

(2.146) 


The sign of A is determined by (2.13). 

Note that in the expansion presented in (2.10) the second sum starts at |£| = N + 1. 
This is due to the fact that we make the additional basis function F(x,y,A) orthogonal to 
e ,kx , thus we use F{x,y,A) — Fiq{x,y,A) in the expansion (2.10). The procedure described 
in (2.14) is second order accurate in the location and jump of the shock. In fact, we can 
state 


Theorem 1: Let u[x ) be a piecewise C°° function with one discontinuity at x,. Let y 
and A be defined in (2.14) then 


|y-z s | =0(jfr) 

\ A - Ml = 0 (ft) • 


(2.15) 


Proof : From (2.3) we get 


(N + 1)uat+i 

e ~i(N+l)x, 

M + ^ + oi 

r i \1 

k(n+i)*J 

[N + 2)u N+2 

e -i{N+ 2)x. 

[[“) + + 0 1 

UAT+l) 3 )] 


= <="•[! + 


By the same token 


|A| — (N + 1)|uat+i 



[u"j Y , M 2 

(iV+l) 2 j (iV + 1) 2 


IMIfi + o(4)]- 
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It should be noted that a better approximation to the shock location x, and its magnitude 
[tt] can be obtained if we add to the basis functions a function of the form 


£ 

\t\>N 


A B 

+ 




(2.16) 


{a {ity j 

and extend (2.11) to \k\ < N + 3. In practice, however, (2.10) is enough to get a nonoscil- 
latory scheme. 

In order to demonstrate that the procedure described in (2.10), (2.12), and (2.14) is 
indeed nonoscillatory we recall the definition of the total variation of a function. 
Definition : The total variation of u over [0, 27 r] - TV[u] - is defined as 


TV[u] = supX)|u(x.) - u(x,_!)| (2.17) 

«=i 

where 0 < xo < zi < • • • < x n < 2n is a partition of [0,27r]. The supremum is taken over 
all partitions. 

It is clear that if u'(x)eL 1 then 

TV[u}= f|«'(0K (2.18) 

Jo 


If we approximate the function u(x) by its finite Fourier series PnU defined in (2.8), then 
it is well-known that the total variation of Pjvu need not approximate that of u. In fact, 
we can state 


Lemma 2: Let the sawtooth function F(x,0, 1) and its N Fourier series Fn(x,0, l) be 
defined by (2.5) and (2.9), then 

TV[F] = 4tt (2.19) 

TV\Fn] = O(logiV). (2.20) 

Proof : Equation (2.19) follows directly from the definition of total variation. As for 
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(2.20) we note that 


2V[F n (z,0,1)] = 

Jo Jo t= _ N 




= / ,r |“ nVL+jh-ndx 

Jo sm t;X 


( 2 . 21 ) 


r2r 

= L i 


2 * .s'm(N + |)x 


sm £x 


| dx + 0(1). 


The first term on the right hand side of (2.21) is the Lebesgue constant. It is known 
([13], p. 67) that it grows like (log N). Hence (2.20). We can therefore conclude that 
TV{Pnu) does not converge to TV[u]. This reflects the existence of large oscillations in 
the neighborhood of the discontinuities. 

The situation is different for v ^ defined in (2.10). In fact, we can state 


Theorem 2: Let u(x) be a piecewise C°° periodic function with one point of disconti- 
nuity x s , and a jump of [u]. Let A and y be such that 

|y - X g \ = A x 

| A - [u]| = A 2 . 

Let vn be defined in (2.10), then 


lo£ N 

TV M < TV[u] + Lq — — + L x A X N log N + L 2 A 2 log N 
— «IU‘ < Cq—jp — ^ l°g^ + C' 2 A 2 . 


( 2 . 22 ) 


(2.23) 

(2.24) 


We present the proof in a series of Lemmas in order to clarify the role of each one of the 
terms on the right hand sides of (2.23) and (2.24). 


Lemma 3: Let F^(x, a, 1) and Fn{x,(3, 1) be defined by (2.6)- (2. 9) and A = a— (3 > 0. 
Then if A < L ; 

TV[F n {x, a, 1) - F n {x, (3 , 1)] = 0(AiV log N) (2.25) 

||F w (x,a,l) - F n (x,P,1)\\ l i = O(AlogN). (2.26) 
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Proof : Since Fn(x, a, 1), F?f(x, /?, 1) are trigonometrical polynomials they are C c 
tions. Therefore, 

TV[F N (x,a,l) - F n (x,P,1)] = f \F' N (x,a,l) - F' N (x,(3,l)\dx 

Jo 


= \ T IE [e il{x ~ a) ~ t il[x ~ p) ]\dx 

, f 2 * if • „l a + P\ ■ , a ~0\j 

= 4 J q |2^sin£(x —)smE—^—\da 


e=i 


Upon defining a t = si we can rewrite (2.27) as 

TV[F N (x,a,l)-F N {x,P,l)} = 4 [*" |f>,sin££|d£ 

Jo e = i 

we note that at are positive and monotone in £ , ot-\ — at < 0. Define now 


to get 


k= 0 


X] o*sin££ = -Bt-i(0) 

i=i t=i 


N 


= — °t)Bt-i + o n b n . 

i=i 


Therefore, 

r2r N 


1=1 

Denote now by /i 
from (2.28) and (2.31) 


[ | J2 a t sin ££ | d£ < a N f \B N (Z)\dt + J2[°i- a t-i) [ |#<-i(£)l rf £* 

JO Jo Jo 




TU[Fftr(a;,a,l) -F N (x,/3, 1)] < 8 fia N . 


In order to estimate fi, we first note that 


B*(£) = sin 


[t+ l)£ _ sinf 
2 sin f ' 


1 func- 

(2.27) 

(2.28) 

(2.29) 

(2.30) 

(2.31) 

(2.32) 

(2.33) 

(2.34) 
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Therefore, 


/>, ( «K< r^fu- 

jo Jo sm £ 


/x<. 


But //* is exactly the Lebesgue constant, therefore, 


(2.35) 


li t = 0(£n£). (2.36) 

Since o N < N A, we get 

TV[F N (x,a,l) -F N (x,(3,l)] = 0{AN\ogN). 

To obtain (2.26), we follow the same arguments as above. Similar to (2.27) and (2.28) 
we have 

f \F N (x,a,l) - F N (x,/3,l)\dx <2nA + 4 [ \ Y] 6*cos££|d£ (2.37) 

Jo Jo f_ x 

where cri = < £ < N. a t s are positive and monotone in £, at — cr*_i < 0. If we define 


Bi{£) = J2 cosk£ = cos ^~ l ~ 

jt=i z 

B 0 {i) = 0 


sinf 
sin | ' 


(2.38) 


Then similar to (2.33) we have 

(It 

/ \Fn(x,ol, 1) — Fff(x,l3,i)\dx < 2nA + 8iiO! (2.39) 

Jo 

where /z is defined in (2.32) with Bt replaced by Bt of (2.38). Notice that (2.35) also holds 
for Bt(x) and \cri\ < A. (2.26) now follows from (2. 35), (2. 37), and(2.39). 


Lemma 4: Let F N (x,a,A ) and Fff(x,/3,B) be defined in (2.6)-(2.9). Denote 


\a-/3\ = A l , \A-B\ = A 2 . 


(2.40) 


Then 


TV[F n {x , a, A) - F n {x, (3,B)]< K^N log N + K 2 A 2 log N 
H-F^x, a,.A) — F n (x,(3,B)\\i,i. < CiAilogN + C 2 A 2 


(2.41) 

(2.42) 



for Ki, K 2 , Ci, C 2 independent of N. 


Proof : 

TV[F N (x,a,A) - F N (x,p,B)} < TV[F N (x,a,A) - F N (x,/3,A)] 

+ TV[F N {x,P,A)-F N {x,P,B)}. 

The first term in the right hand side of (2.37) is bounded by (2.25), the second term 
TV[F n (x,P,A)-F n (x,P,B)\ = T7p-B)£^H 

t-~N 

1*0 


(2.43) 


< \A-B\ F’\ V 

Jo 

< K 2 A 2 \ogN. 


(2.44) 


l=-N 

&0 


Similarly we have 


\\F N {x,a,A) - F n {x,P,B)\\ l i < \\F n (x,oc,A) - F N (x,p,A)\\ L i + \\F n (x,P,A) ■ 

-F N (x,p,B)\\ L i = \\F N (x,oc,A) - F N (x,p,A)\\ L i + \A - B\ ||i^(x,/?,l)|| L i. (2.45) 

The first term on the right side of (2.45) is bounded by (2.26). The second term 


\A-B\ i|^(x,^l)|Ui=A 2 ||^(x,/3,l)|U, <C 2 A 2 


(2.42) follows from (2.45), (2.46), so the Lemma is proven. 


Lemma 5: Let S(x, a, A) and S N (x, a, .A) be defined by 


(2.46) 


oo ii(x-a) 

S(x,a,A) = A 

t =— oo 


{ay 


N e it(x-a) 

S N {x,a,A) = A T, • 


(2.47) 


Then 


TV[S{x,a,A) - 5jv(z,a,A)] < K z 


log N 
N 


(2.48) 
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(2.49) 


for K 3 , Ki independent of N. 

Proof : It is clear that 

fir 

TV{S{x,a,A) — S N (x,a,A)) = / \F(x,a,A) — F N (x,a,A)\dx 

Jo 

the estimates (2.48)-(2.49) follow from [13, p. 185]. 

We are ready now to prove Theorem 2. 

Proof of Theorem 2: First we prove (2.23). In view of (2.3), we can write 

u = F(x, x„ [u]) + S ( x , x„ [u']) + g(x) (2.50) 

and therefore 

P N u = F n (x, x„ [u]) + S N (x, x„ [u']) + P N g{x) (2.51) 

we can also rewrite (2.10) as 

uw(z) = P N u + F(x, y, A) - F n (x, y, A) (2.52) 

hence 

Uiv(x) = FV(x,x s ,[u]) + SW(z, z s , [u']) + P^y(x) + F(x,y,A) - F N {x,y,A) 
or 

Ms) = [-M*. *«. N) - Fn{x, y,A)} + [F(x, y, H) + S (x, y, [«']) + g{x - y + x # )l 
+ [S N (x, x„ [«']) - S (x, y, [«'])] + [P^y(x) - g (x - y + x,)] 

+ [F(x f y,A) -F(x,y, [u])]. 

(2.53) 

The second term on the right hand side is just the original function u shifted 

F(x, y, [«]) + ^(x, y, [«']) + y(x - y + x,) = u(x - y + x s ) (2.54) 
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also from (2.48) 

TV[S n (x, x„ [u']) - S(x,y, [«'])] < TV[Sn(x, x s , [u']) — S N (x, y, [u'])] 

logiV < 2 ' 55 ) 

+ TV[S k {x , y, [«']) - S(l, y, [«'])] < 
and finally since g(x ) is smooth enough 


TV[P N g — g{x — y + x,)] < 


(2.56) 


Therefore from (2.53) and Lemma 4 and Lemma 5 


TVM*)] < TV\u] + + ZqAiJVlogJV + L-A, logW. 


(2.57) 


Next we prove (2.24) following the same argument above. From (2. 50), (2.51), and 
(2.52) 

t>jv(x) - u(x) = [jFjv(x,x„[u]) - FV(x,y,A)] + [S^x,*,, [u']) - S(x,x 4 , [u'])] 

(2.58) 

+ [F{x,y,A ) -F’(x,x i ,[u])3 + [P iV y(x) -y(x)]. 

The first term will be bounded by (2.42), the second term by (2.49), the third term 

( 2 *\F(x,y,A) - F{x,x„[u})\dx < C^Ai + C 2 A 2 . (2.59) 

Jo 

Now since g(x) is smooth enough, we have 

IliVs - dIIl. = O(^). (2.60) 

Therefore from Lemma 4 , Lemma 5 and (2.59)-(2.60) 

|| vn — u||ii < Cq—Jp — f C'lAi logiV + C 2 A 2 , 
and the proof is completed. 


Corollary: The method suggested in (2.15) yields 



and therefore 

TV!v n ] < TVu] + (2.61) 

IK - «||i. < C^. (2.62) 

Thus, the total variation of converges to that of u. converges to u in L 1 norm one 
order higer than P N u for which the rate of its convergence in L 1 is O(^p-). The method, 
therefore, yields a reconstruction technique which is total variation bounded. 

We conclude this section by pointing out that a similar result for collocation method 
and/or for Chebyshev expansions can be developed, along the same lines. Computationally, 
we observe similar results for Galerkin and collocation methods (see Section 4). In practice, 
collocation is used more often than Galerkin, especially when solving a nonlinear PDE 
(Section 3). 


3. NON-OS CILLATORY SPECTRAL SCHEMES 

In this section, we apply the techniques discussed in Section 2 to solve the PDE (1.4): 

u t + /(«)* = 0 (3.1a) 

u(x,0) = u°(x). (3.16) 


If the cell average of u is defined by 

u[x, t) = (2-2) 

then (3.1) can be rewritten as 


^u(x,0 + ^[/(u(x+^. ())-/(< 

II 

O 

(3.3a) 

u(x, 0) = u°(a:). 


(3.36) 


Hence a semi-discrete conservative scheme 



= Mi = -s** i 



(3.4) 
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will be of high order if the numerical flux fj+i approximates f(u(xj + t)) to high order. 

This is the approach used in the MUSCL type semi-discrete finite difference TVD and 
ENO schemes [10], [4]. Notice that (3.4) is a scheme for the cell averages Uj. However, 
in evaluating f- + i, which should approximate f(u(xj + we also need accurate 

point values u ; - + i = u(xj + ^r, t). For finite difference schemes, the reconstruction from 
cell averages to point values is a major issue and causes difficulties, especially in several 
space dimensions [4], [5], For spectral methods, this is very simple because u is just the 
convolution of u with the characteristic function of (Xj_i,x J+ i). To be precise, if 

N 

u ( s ) = a t e%tx (3.5) 

e=-N 

(we have suppressed the time variable t ) , then 

N 

u{x) = ^ ate xtx (3.6) 

<=— JV 

with 

ginf 

a t — a t a t , o t = — j eI— for 0 < Kl < N, a 0 = 1. (3.7) 

2 

Notice that for collocation or Galerkin with Ax = we have |^-| < f for |£| < N, hence 
| < a < 1. The division or multiplication by at thus causes no stability difficulty. We 
point out that at resembles the Lanczos filter [8], which in our notations is 8in ^ 1 ^ , and 
approaches zero when |£| —> N. 

The easy transform between u and u is also valid in several space dimensions and for 
other spectral expansions (e.g., Ghebyshev expansions). We omit the details. 

We now state our scheme as (3.4) with 

f j+ L = f{v N {x j+ L,t)) (3.8) 

where v N is defined by (2.10). We obtain the Fourier coefficients a t of u from {%} by 
collocation, and obtain at of u needed in (2.10) by (3.7). The main difference between the 
conventional spectral method and the current approach is that we use the non-oscillatory 
reconstruction v n instead of the oscillatory P^u in (3.8). 
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The scheme, as it stands, can only treat a solution of not more than one discontinuity. 
However, it can be easily generalized. 

We remark that if u is smooth, (2.10) keeps spectral accuracy because A determined 
by (2.14) will be spectrally small. 

To discretize (3.4) in time we use the high order TVD Runge-Kutta methods in [12]: 


uW = + A*Ai£(uW)], t = 1, • • • ,r 

k =° (3.9) 

u<°) = u fl ) ui r ) = u n+1 . 


In Section 4, we use a third order scheme r = 3 with a 10 = /?io = 1, a 2 o = f, /? 2 o = 


o, «21 = /?21 = <*30 = /?30 = 031 = /?31 = 0, a 32 = /?32 = §• We use a small At so 

that the temporal error can be neglected. These methods are TVD (or TVB) if the Euler 
forward version of (3.4) is TVD (or TVB). In light of (2.61), we expect the total variation 
of (3.4)-(3.8)-(3.9) to grow at most at the rate of 0(inN). In practice we observe stable 
results (Section 4). 

As in the finite difference case [10] [11], we may also apply limiters to obtain provable 
TVB schemes while still keeping spectral accuracy. Let 


Uj — u j+ ± u j i Uj+i — u j+ 1 u j+% (3.10) 

where Uy + 1 = u^(x J+ i ,t) in (3.8). We limit the increments by 

-(mod) _ A + Uy, A_ tty) , Uy +1 ^ = m(u J+1 , A + Uy, A_U ; ) (3.11) 

where m is the minmod function with TVB correction: 

f a u ifjax[<MAx 2 

m(ai, ■ • • , a k ) = < s • mini^* |a,|, if |ai[ > MAx 2 and sign(di) = sVt (3.12) 

( 0, otherwise 

with M = |M 2 or M = Mj = |(3 + 10Af2)M 2 ■ Az i + | A ^!| + | A _ 5 .| • Here M 2 is the maximum 
of |u°J in some region around the smooth critical points of u°(x). See [ll], [2]. 



The flux (3.8) is modified to 



, . ~{mod) _ z{mod) 

h[Uj + u) \u J+1 -u j+l ) 


where h is any monotone flux [3]. We then have the following Lemma. 


(3.13) 


Lemma 6: Scheme (S.9)-(S.1S) is TVB and formally spectrally accurate in space (i.e., 
the spatial local truncation error in smooth region is spectrally small), if the filtering (l.S) 
is used. 


Proof : The proof for TVB can be found in [10], [11]. By [1], the local truncation error 
is spectrally small in smooth regions if the limiter (3.11) returns the first argument. The 
proof that (3.11) always returns the first argument in smooth regions, including at critical 
points, can be found in [2]. 

We remark that the scheme (3.4)-(3.8)-(3.9), with or without the TVB limiting (3.11), 
yields almost identical results in our numerical examples (Section 4). This indicates the 
good stability property of the scheme (3.4)-(3.8)-(3.9). We also remark that (3.13) yields 
a TVB scheme regardless of the underlying method (3.4). However, accuracy in smooth 
regions may be lost if the underlying method (3.4) is globally oscillatory, because the 
limiters (3.11) may be enacted in smooth regions to counter-balance these spurious oscil- 
lations. Numerical examples in Section 4 verify these remarks. In [9], McDonald also used 
some limiters to obtain a TVD spectral scheme. However, the accuracy in smooth regions 
is questionable due to the above remarks. 


4. NUMERICAL RESULTS 

We use several numerical examples to illustrate the methods introduced in the previous 
sections. 


Example 1: We use the approximation (2.10)-(2.12)-(2.14) on the following function 


f sin f , 0 < x < 0.9 

( — sin|, 0.9 < x < 2ir. 


(4.1) 
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Notice that [u^] 7 ^ 0 for all k > 0. Both Galerkin and collocation methods are tested. 
Exponential filters (1.2) with m = 4, k 0 = 0 is used. 

In Table 1, we list the errors of the shock location and shock strength determined by 
(2.14). Notice that the second order accuracy (2.15) is verified. 

Figure 1 displays the numerical solution of Galerkin approximation (2.10)-(2.12)-(2.14) 
with N = 64. Figure 2 is the error of the approximation on a logarithm scale. We have 
found the same kind of results for collocation approximation. In Table 2 , we list the L 1 
error and numerical order in smooth region (in this case, we define the smooth region to 
be 0.8 away from discontinuity). We can see Galerkin and collocation have the same order 
of accuracy. There is no O(l) error near the discontinuity, overall we achieve O(^fr^) for 
L 1 convergence, verifying (2.62). For comparison, we refer the reader to [7]. 

Example 2: We apply (2.10)-(2.12)-(2.14) on a discontinuous function which is the 
steady state solution of an astrophysics problem [ 6 ]. Figure 3 is vn in ( 2 . 10 ) with N = 32. 
For comparison, Figure 4 is the usual Galerkin approximation P N u with N = 32. The 
improvement is apparent. 


Example 3: We solve the Burgers’ equation 



(4.2) 


u(x, 0) = 0.3 + 0.7 sinx 
using scheme (3.4)-(3.8)-(3.9) and (3.4)-(3.9)-(3.13). We find the shock location and 
strength with (2.14). In practice, we do not use the last part of the Fourier modes of 
the numerical solutions in (2.14). In our computation, we find that the coefficients of 
modes in the range of y/N ~ N* give us the best results in shock location and strength. 
It can be proven that in this range of modes (2.14) will not fail in the presence of possible 
transition points in the numerical solutions. The errors of (3.4)-(3.8)-(3.9) in smooth re- 
gions (1.6 away from shock when it appears), at t = 0.8 (before shock), t = 1.42 (when the 
shock just develops), and t = 2.00 (after shock) are listed in Table 3. The numerical solu- 
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tions are displayed in Figures 5-6. The error at t = 2.00, in logarithm scale, is displayed 
in Figure 7. 

We seem to observe higher than algebraic order in smooth regions both before and after 
the shock develops. This might be the first time super-algebraic accuracy is observed in 
a shock capturing spectral scheme solving a nonlinear PDE with shocks. The usual 0(1) 
Gibbs oscillation near the shock is also absent in all of our calculations. We also notice that 
the TVB limiter (3.11) does not change the numerical results significantly in the smooth 
region (see Table 4). Actually, we observe the same order of accuracy in the smooth region, 
comparing Table 4 with Table 3. This indicates that the scheme (3.4)-(3.8)-(3.9) is by itself 
very stable. 

Finally, we run the usual spectral scheme (i.e., with v jy in (3.8) replaced by P^u) with 
the TVB limiter (3.11). The errors in smooth regions (1.6 away from shock) are listed in 
Table 5 (compare with Table 4). Clearly, we only get first order in smooth regions after 
the shock develops. This indicates that TVB limiting can make a scheme stable but may 
not preserve the accuracy. 

Example 4: 2-D Steady State. We solve a 2-dimensional scalar conservation law 
’ u t + ( y) x + u„ = 0 ( x , y)e[ 0, 27 t] x [-1, 1] 

■ u(x, 0, t) = sin 2 (4.3) 

. u(0, y, t) = u( 2 tt, y, t) ye[-l, 1], t > 0. 

We know that (4.3) has a steady state solution Uoo(z,y). u^x, y) actually will be the 
solution to (4.2) if we replace y by t and set u(z,0) = sinx in (4.2). 

As mentioned in Section 3, (3.4)-(3.8)-(3.9) can be extended to 2-dimensional cases, we 
can use either Fourier or Chebyshev method in each of the spatial directions. To solve for 
the steady state of (4.3), we use Fourier method in x-direction and Chebyshev method in 
y-direction. The criteria we set for the steady state is that the relative L 1 residue between 
two consecutive time stage to be less than 10 -6 , i.e., 



Figure 8 displays the profile of the steady state at y = 0.38 and y = 1.00. The solid line is 
the exact solution and the plus signs are the numerical one. 32 points in the x-direction 
and 8 points in the y-direction are used. Figure 8 is the contour plot for the numerical 
steady state solution. 
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Table 1. Errors of shock location and strength, Example 1. 


N 

Galerkin 

Collocation 

Location 

Strength 

Location 

Strength 

Error 

Order 

Error 

Order 

Error 

Order 

Error 

Order 

8 

IggjM 


0.12(-l) 


igigjM 


0.20(-l) 


16 

figfjs 1 

2.6 

0.22 (-2) 

2.4 

llSSpi 

0.8 

0.12(-1) 

0.7 

32 

1 

2.3 

0.48 (-3) 

2.2 


3.8 

0.38(-2) 

1.7 

64 

jtfug 1 

2.1 

0.11 (-3) 

2.1 


2.2 

0.11 (-2) 

1.8 

128 


2.1 

0.27(-4) 

2.1 

BBieii 

2.0 

0.28 (-3) 

1.9 


Table 2. L 1 Error in Region I = {rce[0, 27r] , \x — x 3 \ > 0.8} and Region II = [0, 27r], Example 1. 


N 

Galerkin 

Collocation 

Region I 

Region II 

Region I 

Region II 

Error 

Order 

Error 

Order 

Error 

Order 

Error 

Order 

8 

0.32(-l) 


0.14(0) 


0.23(-l) 


0.3l(-l) 


16 

0.32 (-2) 

3.30 

0.75 (-2) 

4.27 

-0.21 (-2) 

3.46 

0.61 (-2) 

2.34 

32 


3.75 

0.17(-2) 

2.11 

0.23 (-3) 

3.20 

0.17(-2) 

1.79 

64 


5.55 

0.39(-3) 

2.14 

0.54(-5) 

5.40 

0.49(-3) 

1.86 

128 


8.67 

0.96(-4) 

2.04 

0.12(-7) 

8.82 

0.1 3 (-3) 

1.92 


Table 3. Errors in smooth region for (4.2). At t = 0.8, the smooth region is [0, 2 tt], At 
t = 1.42, 2.0, the smooth region is 1.6 away from the shock. 



t = 0 

.8 

t = 1 

42 

t = 2 

.0 

N 

L l Error 

Order 

L 1 Error 

Order 

L 1 Error 

Order 

16 








0.94(-2) 


0.39(-2) 


0.44 (-2) 


32 


3.57 


1.66 


1.40 


0.79(-3) 


0.13 (-2) 


0.16(-2) 


64 


5.00 


5.17 


5.28 


0.25 (-4) 


0.35(-4) 


0.42(-4) 


128 


7.58 


7.71 


6.58 


0.13 (-6) 


0.16(-6) 


0.44(-6) 



22 















































Table 4. Errors in smooth regions for (4.2) of new spectral scheme with TVB limiting 


(3.11). For both t = 1.42 and 2.0, the L 1 errors are taken in the region of 1.6 away from 
the shock. 



t = 1. 

42 

t = 2 

.0 

N 

L 1 Error 

Order 

L 1 Error 

Order 

16 






0.64(-2) 


0.63 (-2) 


32 


1.90 


1.67 


0.17(-2) 


0.20(-2) 


64 


5.55 


5.29 


0.36(-4) 


0.50(-4) 


128 


7.75 


7.11 


0.17(-6) 


0.36(-6) 



Table 5. Errors in smooth regions for (4.2) of the usual spectral scheme with TVB limiting 
(3.11). For both t = 1.42 and 2.0, the L 1 errors are taken in the region of 1.6 away from 
the shock. 



t = 1.42 

t = 2.0 

N 

L 1 Error 

Order 

L 1 Error 

Order 

16 





32 

0.25 (-1) 

0.98 

0.16(-1) 

* 

64 

0.98 (-2) 

1.50 

0.17(-1) 

1.17 

128 

0.34(-2) 

0.84 

0.79(-2) 

0.76 


0.19(-2) 


0.47(-2) 
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Fig. 7: Example 3, Inviscid Burgers’ equation, u(x,0) = 0.3 + 0.7 * sin(x). Errors of 
numerical solutions at time t — 2.00 in the logarithm scale for N = 16,32,64,128. 








Fig. 9: Example 4, Contour plot of the steady solution. 
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